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We discuss real time simulations of high energy nuclear collisions in a classical 
effective theory of QCD at small x. At high transverse momenta, our results match 
the lowest order predictions of pQCD based mini-jet calculations. We discuss 
novel non-perturbative behaviour of the small x modes seen at small transverse 
momenta. 


1 Introduction 


The space-time evolution of nuclei in high energy heavy ion coll^ions and the 
various proposed signatures of the hot and dense^atter formedid depend sen¬ 
sitively on the initial conditions for the evolutionQ. These are the distributions 
of partons in each of the nuclei prior io the collision. In the standard perturba¬ 
tive QCD approach, observables from the collision are computed by convolving 
the known parton distributions of each nucleus with the elementary parton- 
parton scattering cross sections. At the RHIC and LHC colliders, hundreds 
of mini-jets may be formed in the initial collision □. Final state interactior^ 
of these mini-jets are often described in multiple scattering models ^ee Ref. 0 
and references therein) or in classical cascade approaches (see Ref. □ and ref¬ 
erences therein). The possible “quenching” of these mini-jets has also be^ 
studied and proposed as a signature of the formation of a quark gluon plasmau. 
Recently, initial conditions for the energy density and velocity obtained in the 
mini-jet approach have been used in a simple hydrodynamic model to study 
the late time evolution of matter in high energy nuclear collisionsQ. 

In this paper, we describe an ah initio QCD based approach to the study 
of high energy nuclear collisions. This model naturally incorporates coherence 
effects which are important at small x and small transverse momenta and 
simultaneously reproduces the standard mini-jet results at large transverse 
momenta. It also contains a self-consistent space-time picture of the nuclear 
collison. The model is based on an effectivmaction approach to QCD initially 
developed by McLerran and Venugopalan □, and later further developed by 


Ayala, Jalilian-Marian, McLeripj^: 
Kovner, Leonidov and WeigertO’EJ’^^ 


Venugopalanty and by Jalilian-Marian, 
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The above mentioned effective action contains one dimensionful parameter, 
/r, the total color charge squared per unit area. The classical fields correspond¬ 
ing to the saddle point solutions of the effective theory are the non-Abelian 
analogue of the Weizsacker-Williams fields in electrodynamip.^xact analyt¬ 
ical expressions for these fields have been obtained recently Ella. Further, it 
has been shown explicitly that p. obeys renormalization group equations in ra¬ 
pidity y and momentum transfer squared . These reduce to the well Vn nwn 
BFKL and DGLAP equations respectively in the appropriate limits EjEj. 

The above model was applied to the problem of nuclear collisions by 
Kovner, McLerran and Weige^ They formulated the problem as the collision 
of Weizsacker-Williams fields Ej. Further, classical gluon radiation corresond- 
ing to perturbative mod^was studied by these authors and later in greater 
detail by several othersEj’Ej’Ej. In the small x limit, the classical Yang-Mills £&- 
suit agrees with the quantum Bremsstrahlung result of Gunion and BertschcJ. 

The perturbative approach is of course very relevant and useful. However, 
it is essential to consider the full non-perturbative approach for the following 
reasons. Firstly, the classical gluon radiation computed perturbatively is in¬ 
frared singular, as is also the case in mini-jet calculations, and therefore rather 
sensitive to the cut-off. It was argued in Ref. E§EZl that a natural scale for 
the cut-off is kt ~ ag/r. However, to be sure, it is important to perform a 
full calculation. Secondly, the non-perturbative approach is crucial to study 
the space-time evolution of the nuclei and in particular, the possible thermal- 
ization of the system. This has several ramifications for various quark gluon 
plasma signatures. For instance, if thermalizjation does occur, then hydrody¬ 
namic evolution of the system is reasonable Eu. In that event, our approach 
would provide the necessary initial temperature and velocity profiles for such 
an evolutions. 

We discuss in this paper first results from real time simulations of the 
full, non-perturbative evolution of classical non-Abelian Weizsacker-Williams 
fields. Such a simulation is possible since the fields are classical. Similar 
classical simulations have ^en performed in the context of sphaleron induced 
baryon number violation E3 and chirality vic^ting transitions in hot gauge 
theories CJ. In brief, the idea is as follows a. We write down the lattice 
Hamiltonian which describes the evolution of the small x classical gauge fields. 
It is the Kogut-Susskind Hamiltonian in 2-|-1-dimensions coupled to an adjoint 
scalar field. The lattice equations of motion for the fields are determined 
straightforwardly. The initial conditions for the evolution are the Weizsacker- 
Williams fields for the nuclei before the collision. Interestingly, the dependence 
on the sources does not enter through the Hamiltonian but instead through 
the initial conditions. 
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A related approach is that of Mueller, Kovchegov and WallonB, where they 
combined Mueller’s dipole picture of high energy scatteringc3 with the classical 
Yang-Mills picture O to study nucleon-nucleus scattering. For alternative 
approacheSpWe refer the reader to the work of Makhlin and SurdutovichcHand 
of BalitskyS 

Our paper is organized as follows. In the following section we briefly discuss 
the problem of initial conditions for nuclear collisions and the perturbative 
computation of classical gluon production. We also discuss a non-perturbative 
Hamiltonian approach to the solution of the full Yang-Mills equations. In 
section 3, we formulate the problem of solving the Yang-Mills equations on 
the lattice. Starting from the lattice action and assuming boost invariance, 
we write down the lattice Hamiltonian, the Hamilton equations of motion and 
the initial conditions for the dynamical fields and their conjugate momenta on 
the lattice. The initial conditions depend on a single dimensional parameter 
/i, the color charge per unit area per unit rapidity. Numerical results from 
our simulations and comparisons to lattice perturbation theory (LPTh) are 
discussed in section 4. These are performed for a range of values of /x = 0.015- 
0.2, and for lattice sizes from 10 x 10 to 160 x 160 measured in units of the 
lattice spacing. We summarize our results in section 5. 

2 The Weizsacker Williams approach to high energy nuclear colli¬ 
sions 

The model of McLerran and Venugopalan of gluon fields in a nucleus atpmall 
X was applied to nuclear collisions by Kovner, McLerran and WeigertES. We 
shall now discuss their formulation of the problem and their perturbative com¬ 
putation, to second order in the parameter ag/x/fct, of classical gluon radiation 
in nuclear collisions. We will then briefly discuss a non-perturbative Hamil¬ 
tonian approach which suggests how all orders in as^l/kt can be computed 
numerically. The implementation of this approach on the lattice is described 
in section 3. 

In high energy nuclear collisions, hard valence parton modes are highly 
Lorentz contracted, static sources of color charge for the wee Weizsacker- 
Williams modes in the nuclei. The sources are described by the current 

= 5’^+gpl{n)5{x-) + 5’^-gp^_{n)5{x +), (1) 

where /9+,- are the color charge densities of the hard modes. The classical 
field of the two nuclei is obtained by solving = J'^ for the two above 

mentioned light cone sources. 
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Gluon distributions are simply related to the Fourier transform Af{kt) 
of the classical field by {Af{kt)Af{kt))p. The subscript p above refers to the 
averaging over the color charge distributions (performed independently for each 
nucleus) with the Gaussian weight p?. We assume equal A for simplicity. 

We have omitted the rapidity dependence of the charge distributions in the 
equations immediately above. The omission is discussed further in the next 
section. We note her^hat the rapidity dependence of the charge distribution 
is also absent in Ref.Ej (see the discussion below Eq. ^). 

Before the nuclei collide (t < 0), a solution of the equationa^of motion is 
A^ = 0, A® = d{x~)9{—x~^)a\{rt) + 9{x'^)9{—x~)a2{rt) where, 0 


>'i.2(ri) = — 
ig 


Pe 


-^^11 




dr]' P±W ,rt) 


t 


V® Pe 


—ig d 7 ]'^ 

^ J±Vr,T-^-i ^ , 


P±(v',r-t) 


( 2 ) 

Above, r] = rjproj — log(a; Ixproj) i® the rapidity of the nucleus moving along 
the positive light cone with the gluon field a\ and p = —rjproj + log(a;)J(,jj^ /a;“'') 
is the rapidity of the nucleus moving along the negative light cone with the 
gluon field a\. At central rapidities, (or a: <C 1) the source density varies slowly 
as a function of rapidity and a® = Q;®(rt). The above expression suggests that 
for t < 0 the solution is simply the sum of two disconnected pure gauges. Just 
as in the Weizsacker-Williams limit in QED, the transverse components of the 
electric field are highly singular. 

For t > 0 the solution is no longer pure gauge. Working in the Schwinger 
gauge x''^A~ + a;“A+ = 0, the authors of Ref. E-3 found that with the ansatz 
A^ = ±x^a{T,rt), A® = a(^(r, r^), (where r = V‘2x~^x~), the Yang-Mills 
equations could be written in the simpler form 


^drT^dra + [Di, [iJ®,Q;]] = 0, 

-[Di,dra\]+igT[a,dra] = 0 , 

T 

—drTdra\ — igT^[a, [iJ®, a]] — [D^, pA] = 0 . (3) 


Note that the above equations of motion are independent of ij-the gauge fields 
in the forward light cone are therefore only functions of r and rj and are 
explicitly boost invariant. 

The initial conditions for the fields a(r, r*) and at r = 0 are obtained 
by matching the equations of motion at the point x^ = 0 and along the 
boundaries a;+ = 0, > 0 and x~ = 0,a;+ > 0. Because the sources are 

highly singular functions along their respective light cones, so too in general 
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Figure 1: The lattice size dependence of the scalar kinetic energy density, expressed in units 
of for fi = 0.025 (pluses) and fi = 0.05 (diamonds). The solid line is the LPTh prediction. 
The error bars are smaller than the plotting symbols. 


will be the equations of motion. Remarkably, there exist a set of non-singular 
initial conditions for the smooth evolution of the classical fields in the forward 
light cone. In terms of the fields of each of the nuclei before the collision 
{t < 0), they are, ct^|=_o = cul + and q;|t-=o = Gyulassy and 

McLerran have shown Ej that even when al 2 are smeared out in rapidity to 
properly account for singular contact terms in the equations of motion, the 
above boundary conditions remain unchanged. Further, since the equations 
are very singular at r = 0, the only condition on the derivatives of the fields 
that leads to regular solutions are 9Ta|T=o, 9 tQ!^|t=o = 0. 

In Ref. 113 , perturbative solutions (for small p) were found to order by 
expanding the initial conditions and the fields in powers p (or equivalently, in 
powers of asp/h) as a = ^ “!(«)• subscript n 

denotes the nth order in p. For details of the perturbative solution, we referee 
reader to their paper (see also the related paper of Gyulassy and McLerrantj). 

The perturbative Yang-Mills result for the number distribution of classical 
gluons is 


dN 

dy(Pkt ( 27 r)‘^ kf 




( 4 ) 
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where L{kt,X) is an infrared divergent function at the scale A. We first note 
that thi^esult agrees with the quanturn-hremsstrahlung formula of Gunion and 
BertschEd and with several later works ESEZlilS. It was also shown by Gyulassy 
and McLerran that for sources smeared out in rapidity, the resulting expression 
is identical to the one above except ^ X~^{y)x~iv)- The ± superscripts 
refer to the nucleus on the positive or negative light cone respectively. 

The origin of the infrared divergent function L(kt,X) above is from long 
rang^olor correlations that are cut-off eitheiiiy a nuclear form factor (as in 
Refs.E3’EZl),±iv dynamical screening effects E3 ’Ej or in the classical Yang-Mills 
case of Ref.Ej, non-linearities that become large at the scale kt cesfJ,- In the 
classical case then, L{kt,X) = log{kf/X‘^) where A = ag/r. At sufficiently high 
energies, the behaviour of gluon radiation infrared is given by higher order (in 
asy/kt) non-linear terms in the classical effective theory. 

It is unlikely that a simple analytical solution exists for the Yang-Mills 
equations in general. The classical solutions then have to be determined nu¬ 
merically for t > 0. The straightforward procedure would be to discretize the 
Yang-Mills equations but it will be more convenient for our purposes to con¬ 
struct the lattice Hamiltonian and obtain the lattice equations of motion from 
Hamilton’s equations. This will be done in the next section. Before we do 
that, we will discuss here the form of the continuum Hamiltonian. 



Figure 2: Normalized field intensity of a hard (kt = 2.16GeV) mode vs proper time t in 
units of fm (diamonds). Solid line is the LPTh prediction. 
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We start from the QCD action (without dynamical quarks) In the forward 
light cone {t > 0) it is convenient to work with the r, rj, ft co-ordinates where 
T = \/2x+x~ is the proper time, ij = ^ log(a;“''/a;“) is the space-time rapidity 
and ft = (a;, y) are the two transverse Euclidean co-ordinates. In these co¬ 
ordinates, the metric is diagonal with = 1 and g^^ = 

After a little algebra, the Hamiltonian can be written asH 


H = j dgdftT + ^Pr,rFr,r + jF^yFa^y + f Ay + f Ar^ . 

( 5 ) 

Here we have adopted the Schwinger gauge condition which is equivalent to 
requiring A^ = 0. Also, p"^ = ^drAy and p^ = rdrA^ are the conjugate 
momenta. 

Consider the field strength Fy^ in the above Hamiltonian. If we assume 
approximate boost invariance, or Ar{T,ri,Ft) « Ar(T,Ft) ; Ay{T,g,Ft) ~ 
d)(T,ry), we obtain Fy^ = —IIr‘I*“. Here Dr = dr — igAr is the covariant 
derivative. Further, if we express in terms of the defined in Eq. we 
obtain the result that = 0 for r > 0. Finally, since = T^a(r, ?y) ; 
Ar = a(j_(r, 7y), we can perform the integration over the space-time rapidity 
to re-write the Hamiltonian in Eq. as 


H = 


dr try - 


1 / i9a’l 


dr 


\K,K, + l^{Ds[r^a])- 


( 6 ) 


Here the index (3 = {r^fl). The excrete version of the above Hamiltonian is 
the Kogut-Susskind Hamiltonian t3 in 2+1-dimensions coupled to an adjoint 
scalar field. The lattice Hamiltonian will be discussed further in the next 
section. 

We now briefly comment on a key assumption in the above derivation, 
namely, the boost invariance of the fields. This invariance results in Eq. 
thereby allowing us to restrict ourselves to a transverse lattice alone. There has 
been some confusion about the importance of boost invariance as a consequence 
of the initial assumption that the classical currents a^^lta functions on the 
light cone. For remarks clarifying this issue, see ref. E2l’E3. In general, at the 
energies of interest, particle distributions are unlikely to be boost invariant □. 


3 Real time lattice description of nuclear collisions 


In this section, we discuss the problem on nuclear collisions on the lattice. 
Taking the naive continuum limit in the longitudinal directions of the full 4 
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Figure 3: Normalized field intensity of a soft (kt = lOSMeV) mode vs proper time t (in 
units of fm) for fj, = 200MeV (squares), /i = lOOMeV (pluses), and fi = SOMeV (diamonds). 
Solid line, nearly coinciding with the fi = SOMeV curve, is the LPTh prediction. 


dimensional Minkowski Wilson action for the SU(Nc) gauge group in funda¬ 
mental representation, one obtains an expression for the action S', which-is 
discrete in the transverse directions and continuous in the z-t directions E3. 
The equation of motion for a field is obtained by varying S with respect to 
that field. Initial conditions on the lattice are derived from the lattice ana¬ 
logue of the continuum initial conditions. We start from the lattice action S, 
obtain the lattice equations of motion in the four light cone regions and de¬ 
termine non-singular initial conditions by matching at t = 0 the coefficients 
of the most singular terms in the equations of motion. The full derivation is 
described in Ref.Ej. 

The initial condition (at t — 0) for the transverse link field for the SU(2) 
case is 


U = 


Tr(I/i + U 2 ) 

ul + ul 


/ = ([/! + U2){Ul + Ul)-^. 


(7) 


For Nc > 2 the solution for U is not as simple. This condition has the right 
formal continuum limit. This can be seen by writing Ui ^2 as exp(w^ai_ 2 ) and 
U as exp(m_LQ;j_)- One obtains a±_ = ai -I- Q !2 as required. 

The result for the rapidity component of the gauge field (recall that we 
































































had = ±x^a{xt,T)) is 

= 4^ E -I)- h.c.],,. 

n 

- [{U^ - I){Ui - U 2 ) - h.c.],_„,„) . (8) 

As above, in the limit of smooth fields, one obtains as required a = i '^ 2 ]"- 

The lattice Hamiltoniajiis obtained from the action following the standard 
Kogut-Susskind procedures, giving 



+ ^ETr ($, -C/,(9) 

h" j 

where Ei are generators of right covariant derivatives on the group and Uj^n is 
a component of the usual SU(2) matrices corresponding to a link from the site 
j in the direction n. The first two terms correspond to the contributions to 
the Hamiltonian from the chromoelectric and chromomagnetic field strengths 
respectively. In the last equation $ = $“ 17 “ is the adjoint scalar field with its 
conjugate momentum p = Matching the equations of motion at r = 0 

gives us the initial conditions Ei\t=o = 0 , d)jjT-=o = 0 and p_;jr=o = 2 a. 

Lattice equations of motion follow directly from Hl of Eq. H. For any 
dynamical variable v with no explicit time dependence, v = {Hl,v}, where v 
is the derivative with respect to r, and {} denote Poisson brackets. We take Ei, 
Ui, Pj, and <i>j as independent dynamical variables, whose only nonvanishing 
Poisson brackets are 

{pt, $5} = S,,Sab; {Er, U^} = -iSimUia^; {Ef, E^} = 25i^eabcEf 

(no summing of repeated indices). The equations of motion are consistent with 
a set of local constraints (Gauss’ laws). 

4 Mini—jets on the lattice: numerical results and comparison to 
lattice perturbation theory 

In this section, we will compare results from our numerical simulations to 
analytic results from lattice perturbation theory. Before we do that we should 
first understand the ramifications of our results for nuclear collisions. 
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Figure 4: Field intensity over as a function of kt for fj, = 200MeV (squares), ^ = lOOMeV 
(pluses), and p = 50MeV (diamonds). Solid line is the LPTh prediction. The field intensity 
is in arbitrary units and kt is in GeV. 


In section 2, we introduced the scale which is the color chaxge squared 
per unit area. Its magnitude can be determined from the relationllS 

where q, g stand for the nucleon quark and gluon structure functions at the 
resolution scale Q of the physical process of interest. Also, above xq = Qj\fs. 
Using the HERA structure function data, Gyulassy and McLerran estimated 
that /X < 1 GeV for LHC energies and g < 0..5 GeV at RHIC. Thus the 
regime where the classical Yang-Mills picture can be applied is rather limited. 
However, a window of applicability does exist and depending on what higher 
order calculations will tell us, this window may be larger or smaller than the 
naive classical extimates. _ 

Gyulassy and McLerran Ej have shown that the classical Yang-Mills for¬ 
mula in Eg, is at small x (approximately) the same as the perturbative QCD 
predictionB for the process AA —s- g 


da- _ asNc f ,2 f{xi,qf)f{x2, {h - qtf) 

dycPkt ^ J 


( 11 ) 
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where f{x,Q‘^) = and xi « X 2 = kt/y/s. The two formulae are 

equivalent if we divide the above formula by ttR^, approximate the integral 
above by factoring out / above at the scale and taking the normalization 
factor Kn 5. 

At large transverse momenta, the field intensity measured on the lattice, 
|A(fct,r)p, is simply related to the Yang-Mills distribution function. Then, 
from the above discussion, the field intensity of the hard modes on the lattice 
can also be simply related to the mini-jet cross section for the process AA —> g. 
However, at smaller transverse momenta, there is no simple relation between 
the field intensity and the classical gluon distribution function (and thereby 
the cross section by the above arguments). We have to look for more general 
quantities which, conversely, in the limit of large kt , will ave us the AA g 
cross section. This point will be discussed further in Ref. E3. 



Figure 5: Time history of the energy density in units of for /r = 200MeV (squares), 
fj, = lOOMeV (pluses), and fj, = 50MeV (diamonds). Error bars are smaller than the plotting 
symbols. Proper time r is in fm. 


We now turn to numerical results from our simulations. These were per¬ 
formed for a variety of lattice sizes L= 20-160 and the color charge density 
/t = 0.025-0.3 in units of the lattice spacing a. To convert lattice results to 
physical units, we take « ttR^ for A=200 nuclei. This then also determines 
fjL for a fixed lattice size. The relation of lattice time to continuum time is 
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given by the relation tc = a • tl In Fig. |^, we plot the Gaussian averaged 
initial kinetic energy {Ek) on the lattice as a function of the lattice size L and 
compare it with the lattice perturbation theory expression 


= N,{Nl - 1 ) 


( 12 ) 


E 


E 


sin(l„) sin(Z„/) 

A ^(0 


iME 


sin^(T)sin^( 2 ' 

A2(Z) 


where A(^) = 2 ~ cos(27rZ„/L)], n = x,y, and 1 — L/2 < lj;^y < Lj^. The 

comparison is made for values /i = 0.025,0.05 of the color charge density. For 
small values of L, there is very good agreement between the two but for the 
largest value L — 160, they begin to deviate. The strong coupling parameter 
on the lattice is oc and for ^ 1, we can expect to see deviations 

from lattice perturbation theory. 

In Fig. we plot the ratio of the field intensity of a particular mode of the 
transverse gauge field as a function of proper time r normalized to its value 
at r = 0. The diamonds are results from a lattice simulation with L = 160 
and g = 0.025 and the mode considered is {kx,ky) = (7r/4,0). Note that 
kx,y = 2TTlx,y/L and for this case, lx = 20, ly = 0. The solid line in the fig¬ 
ure is the square of the Bessel function Jq{ujt) where uj = A{1). The time 

dependence of the high transverse momentum modes should agree with the 
continuum perturbative result which predicts a time dependence proportional 
to Jq (ojt) for the field intensity of the transverse gauge fields. The continuum 
dispersion relation uj = |fcj_| is however modified into the above mentioned lat¬ 
tice dispersion relation. We see from the figure that the anticipated agreement 
between the lattice results and perturbation theory is quite good. 

In Fig. 1^, we plot the same quantity as in Fig. 2, but now for three different 
values of g, and for the first non-zero momentum mode {kx, ky) = (1,0). The 
lattice size L = 160 is the same as previously. For the smallest value of 
fjL = 0.025, there is again an agreement with the Bessel behaviour predicted 
by perturbation theory. However at the larger values of /r = 0.05,0.1, one 
sees significant deviations away from the Bessel behaviour. Indeed the modes 
appear to saturate at larger values of r. It is not clear that this saturation has 
much significance since the energy shows the expected 1/r behaviour at late 
times. The late time behaviour of long wavelengths may be understood better 
by looking at gauge invariant quantitiescj. 

In Fig. the field intensity of the transverse gauge field normalized by 
at r = 0 is plotted as a function of the transverse momentum in physical units. 
The lattice results for the different values of g described in the caption are 
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compared to the lattice perturbation theory result below for the field intensity 
(the prime denotes that a satisfies the Coulomb gauge condition): 

, g^iVc(iVe^ - 1 )m" ^ A(2/^ - 0A(0 - [A(r - 0 - 
' 4A(0 ^ A^l')A^l'-l) ■ '' ’ 

The LPTh result (which would be the mini-jet distribution in the continuum) 
agrees very well with the lattice result for small /i up to very small values of 
kt- However, strong coupling effects grow with increasing fj, (the lattice size 
L is fixed here) and we see deviations from the perturbative predictions at 
larger values of kt- This trend is enhanced further at larger values of /r than 
those shown here. The non-perturbative effects due to the non-linearities in 
the Yang-Mills equations seem to temper the 1/kf behaviour predicted by 
perturbation theory. Whether this reflects the presence of a time-dependent 
mass in the theory needs further investigation. 

Finally, we plot the time dependence of the energy density for different 
values of /r. At late times, from general considerations we expect that E oc 
l/r and that is indeed what we see. Perturbation theory predicts the energy 
densities will scale as Clearly our results do not show this scaling which 
suggests that non-perturbative effects are important. It appears though that 
the trend as we go to smaller values of /i is to approach the perturbative scaling 
behaviour. 

5 Summary and outlook 

In this paper, we have described results from real time lattice simulations of the 
full classical Yang-Mills equations with initial conditions given by the classical 
fields of each nucleus before the collision. At large transverse momenta, our 
simulations agree very well with lattice perturbation theory as they should. At 
small transverse momenta, we show that there are significant deviations from 
the LPTh predictions. The late time behaviour of these modes is particularly 
interesting. The spatial dependence of equal time correlators provide further 
insight into these modes. These will be discussed in a forthcoming paper E3. 
While this work was being prepared, we received a preprint studying lattice 
simulations of nuclei in a different approach l3. 
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